%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Discrimination in Multi-Phase Systems: Evidence from Child Protection 

%Created: 2/8/2023 
%Updated: 6/6/2023 

%This program estimates the Hierarchical MTE Model in the paper, following 
%closely the approach in Arnold Dobbie Hull (AER, 2022) 

%Note that we have removed the file directory names from this program for confidentiality reasons.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% SET UP %%%%%%%%
clear all;
path=; 
addpath('subroutines/');

%Options
time='';          
R_init=100;          
R_final=500;         

options_init=struct('Display','iter','MaxFunEvals',100000,'MaxIter',5000,'GradObj','off','TolFun',1e-3);
options_final=struct('Display','iter','MaxFunEvals',10000000,'MaxIter',50000,'GradObj','off','TolFun',1e-5);

%% indexing(row,col) : means everything 
for nosigma=0:0 
    for nobeta=1:1 
        
         %%%%%%% BRING IN DATA %%%%%%%% 
         b_ests=csvread(strcat(path,'data/',time,'b_ests_qje.csv'),1); %% bring in estimates
         b_ests=[b_ests(:,1);b_ests(:,2)]; %% rowbinding b_ests(:,1) on top row and (:,2) on bottom row
         L=size(b_ests,1)/2; 
         w_ests=csvread(strcat(path,'data/',time,'w_ests_qje.csv'),1);
         w_ests=[w_ests(:,1);w_ests(:,2)];
         ests=[b_ests,w_ests]; %% column appending
         
         b_var=csvread(strcat(path,'data/',time,'b_var_qje.csv'),1);
         b_Omega=[b_var(:,1);b_var(:,2)];
         b_Omega=diag(b_Omega); %% only keep the diagonals
         w_var=csvread(strcat(path,'data/',time,'w_var_qje.csv'),1);
         w_Omega=[w_var(:,1);w_var(:,2)];
         w_Omega=diag(w_Omega);

         %% inverting parts of matrices (after rowbinding them)
         b_W_fs=inv(b_Omega(1:L,1:L)); w_W_fs=inv(w_Omega(1:L,1:L));
         b_W_ss=inv(b_Omega(L+1:2*L,L+1:2*L)); w_W_ss=inv(w_Omega(L+1:2*L,L+1:2*L));
         W={b_W_fs,b_W_ss,w_W_fs,w_W_ss}; %% creating new array by appending

        %%%%%% ESTIMATE %%%%%%
        rng(42); %% set the seed for replication

        %% creating random draws from normal with mean 0 and sd 1 and creates a [row size,column size] matrix
        b_param_sims_init=normrnd(0,1,[2*L,R_init]); w_param_sims_init=normrnd(0,1,[2*L,R_init]);
        b_param_sims_final=normrnd(0,1,[2*L,R_final]); w_param_sims_final=normrnd(0,1,[2*L,R_final]);
        b_error_sims_init=normrnd(0,1,[2*L,R_init]); w_error_sims_init=normrnd(0,1,[2*L,R_init]);
        b_error_sims_final=normrnd(0,1,[2*L,R_final]); w_error_sims_final=normrnd(0,1,[2*L,R_final]);

        sims={b_param_sims_init,b_param_sims_final,b_error_sims_init,b_error_sims_final,w_param_sims_init,w_param_sims_final,w_error_sims_init,w_error_sims_final};

        theta_0=[-0.5,0.5,0,-0.5,1,0,-0.5,0.5,0,-0.5,1,0,0,0]'; 

        if nobeta==1 && nosigma==0 
            theta_0=[theta_0(1:2);theta_0(4:8);theta_0(10:14)]; 
        end 
        if nobeta==0 && nosigma==1
            theta_0=[theta_0(1:3);theta_0(5:9);theta_0(11:13)]; 
        end 
        if nobeta==1 && nosigma==1 
            theta_0=[theta_0(1:3);theta_0(5:9);theta_0(11:13)]; 
        end   

        [theta_hat_final,theta_var,theta_hat_init,Q]=get_estimates(ests,sims,b_Omega,w_Omega,W,theta_0,options_init,options_final,nobeta,nosigma);
        theta_se=sqrt(diag(theta_var)); 
        [theta_hat_init,theta_hat_final,theta_se] 
        
        %% saving results
        save(strcat(path,'data/',time,'both_nobeta',num2str(nobeta),'_nosigma',num2str(nosigma),'_ests_qje.mat'),'theta_hat_final','theta_hat_init');
        save(strcat(path,'data/',time,'both_nobeta',num2str(nobeta),'_nosigma',num2str(nosigma),'_var_qje.mat'),'theta_var');

        %end
        
    end
end
